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Abstract 

A two-temperature linear spin model is presented that allows an easily understandable intro- 
duction to non-equilibrium statistical physics. The model is one that includes the concepts that 
are typical of more realistic non-equilibrium models but that allows straightforward steady state 
solutions and, for small systems, development of the full time dependence for configuration proba- 
bilities. The model is easily accessible to upper-level undergraduate students, and also provides a 
good check for computer models of larger systems. 



1 



I. INTRODUCTION 



For over a century, the statistical mechanics of systems in thermal equilibrium has been 
a well-established part of the physics core curriculum. Yet, for natural systems, thermal 
equilibrium is the exception rather than the rule. If we look around us, we see a world 
overwhelmingly far from equilibrium: from the intricate dynamics of living cells to more 
complex biological organisms; from ripples on a pond to global weather patterns. 

The study of systems far from thermal equilibrium is both challenging and rewarding. The 
reward comes from the chance to better understand the vast collection of non-equilibrium 
real-life systems. The challenge is obvious: we step outside the comfort of the familiar 
Gibbs 1 framework for equilibrium systems, so we have to search for new tools to achieve 
fundamental understanding and theoretical classification of non-equilibrium behavior. Over 
the last three decades, an increasing number of condensed matter theorists are devoting 
their efforts to understanding complex collective behavior of far-from-equilibrium systems 
using methods that range from easily accessible computer simulations to sophisticated field 
theoretic techniques. 

Non-equilibrium statistical physics is a relatively new field, not yet part of the physics cur- 
riculum at the undergraduate level. Nonetheless, the methods and results of non-equilibrium 
statistical physics are being employed in fields as diverse as molecular biology, computer 
science, economics and politics. Thus it seems suitable for undergraduate students of quan- 
titative science to be introduced to the field early in their academic journey. 

This paper presents an introduction to the field of non-equilibrium statistical physics 
via a theoretical model which, while retaining the essence of the difficulties of far-from- 
equilibrium systems, is as simple as possible. The "two temperature kinetic Ising model" 
has two appealing features: multi-temperature systems are fairly common, for instance a 
water tank with an immersion heater, nuclear magnetic resonance in an external magnetic 
field, or a lattice of nuclei in a solid prepared at a certain spin temperature^; as well, the 
model can be solved analytically for small system sizes using only basic linear algebra, and 
can also be modeled via Monte Carlo simulations without too much difficulty. For these 
reasons, this model has great pedagogical value. Undergraduate students in a traditional 
course in statistical physics can extend principles learned there to gain first exposure to 
far-from-equilibrium systems. 
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This paper is structured as follows: We first (Section (XT]) give an overview of non- 
equilibrium statistical physics and compare it with its equilibrium counterpart. Next, (Sec- 
tion [III]) we introduce the "two temperature kinetic Ising model" as a simple example of 
a non-equilibrium system. Using the microscopic approach, we solve exactly the master 
equation for some small system sizes, to obtain the probability distribution (Section HVT) . 
We also comment on the probability currents (Section [V}, a fundamental feature of far 
from equilibrium systems in the steady state. By direct calculation, we also find the time 
dependence of the probability distribution for a small system (Section |VT|) . 

II. IN AND OUT OF EQUILIBRIUM: COMPARISON AND CONTRAST 

Statistical mechanics is the discipline which bridges the gap between the microscopic 
world of molecules, atoms and electrons and the macroscopic world of thermodynamics 
and bulk properties of materials. It enables us to predict macroscopic behavior of a system, 
knowing the quantitative properties of molecular interactions. The foundation of equilibrium 
statistical mechanics rests upon Boltzmann's fundamental hypothesis 3 . 

If an isolated macroscopic system is ergodio 4 , it will reach thermal equilibrium 
after a sufficiently long time. Then every configuration (or microstate) available 
to the system can be found with equal probability. 

Of course, completely isolated systems are rather unrealistic, so typically we consider systems 
in contact with a very large (infinite) heat reservoir. After a sufficiently long time equilibrium 
is reached, meaning that the average net energy flux between the system and the thermal 
bath vanishes and both have reached the same temperature, T. Under these conditions, the 
probability for finding the system in a microstate a is given by the canonical distribution: 



where j3 = l/(ksT) is related to the inverse temperature via Boltzmann's constant ks, 
and Z is the factor that ensures the normalization of probabilities, known as the partition 
function. Thus, once we have specified a labeling of the microscopic configurations, {a}, and 
we have determined the microscopic Hamiltonian H(a) — i.e. the internal energy of each 
configuration — so that the Boltzmann factor, e~^ n , is known, we can calculate, at least in 
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principle, the partition function of the equilibrium system as well as average values of time- 
independent observables. Of course, we can run into technical difficulties - in particular, 
some of the configurational sums may not be obtainable exactly, etc. - but the fundamental 
framework for a solution is fully established. 

The jump from the idealization of thermal equilibrium to the full, possibly turbulent, 
dynamics of the real world is too great for our limited present knowledge of complex non- 
equilibrium systems, and we are therefore well advised to focus on the simplest extension 
of equilibrium systems, namely, non- equilibrium steady states. This category of systems 
is characterized by time-independent macroscopic behavior which results from applying a 
uniform driving force. We may model the external driver as a second temperature bath, at 
a different temperature than the first, thus feeding energy into the system or removing it. 
The long-time behavior of such a two-bath system exhibits constant energy flow through 
the system. A simple example of such a system is a curent-carrying electrical resistor in a 
steady state, gaining energy from a current source and losing it as heat to the environment. 
The resistor is not in equilibrium because there is a non-zero energy flux flowing through 
it; yet, after a sufficiently long time, it reaches a "steady state", with time-independent 
macroscopics. Understanding such a system in the steady state involves finding the associ- 
ated stationary probability distribution of microstates, which is the long-time limit of the 
time-dependent microstate probability distribution. 

A starting point for the study of these systems is the master equation, expressing the 
conservation of configurational probabilities. We consider a continuous-time dynamics, with 
a finite and discrete configuration space a. The time-dependent probability Pio, i) for 
finding the system in configuration o at time t changes only due to transfer of probability 
into a from other configurations, or from a into others, in such a way that P(a,t) = 1 
at all times. The evolution of probability P(cr,t) is dictated by a set of transition rates 
c[a — > er'] that describe the evolution of the system from one configuration a to a different 
configuration a' , per unit time. For example, given a spin system, one configuration leads 
into another via a spin flip. We may write a balance (continuity) equation: its right hand 
side consists of two sums: the first is a "gain" term, summing over all configurations from 
which configuration a could possibly result while the second is a "loss" sum, accounting for 
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all configurations into which a can evolve: 

= E < c v - ff ] *v> *) - c [* - ^i ^ 0} (2) 

a' 

Here, c [cr — > a 1 denotes the transition rate from configuration cr into another configu- 
ration cr'. These rates must be given as part of defining a specific model. Our task 
in solving the steady-state problem is to find the stationary solution of this equation, 
P*(cr) = linit-,.00 P((T, t), for which the left hand side of Eqn. (j2J) vanishes: 

= d -^fl = J2a> {C W - a] P*(a) - C [cr -> cr'] P»} (3) 

Of course, the steady state distribution P*(cr) will depend on the transition rates. Under 
rather general conditions on the c's, this solution will be unique, and thus independent of 
initial conditions. 

For a system in thermal equilibrium with a single heat bath, we know its steady state 
distribution must be given by P eg (cr) from Eq. (JTJ). For the equilibrium case, we must thus 
choose the rate terms c [cr — > cr'] to be consistent with the canonical distribution result. This 
requirement forces the rates to satisfy the "detailed balance" condition, 

C [cr' -> a] _ P e(? (cr) 



c[a —> <j'] PeqW) 



(4) 



Since P eq (cr) oc exp(—/3H), we therefore choose our rates such that 

= exp(/3AP0 (5) 

c[a — > cr J 

where 

A# = tf(cr') - #(<t) (6) 

For systems in equilibrium, the "detailed balance" condition is an intrinsic property, and 
is related to their microscopic reversibility 5 . This condition assures the invariance of the long 
time limit of the probability distribution. This constraint on the transition rates is necessary 
when modeling (via Monte Carlo simulations, for instance) systems in thermal equilibrium. 
A key feature of a system far from thermal equilibrium is the violation of detailed balance: 
its steady state distribution P*(cr) does not satisfy Eqn. (p|J. For a non-equilibrium system 
in its steady state, the question becomes how one generalizes the detailed balance condition 
such that, when the "drive" is turned off, the equilibrium solution is recovered. 



A more intuitive way to discuss the detailed balance condition is to describe it in terms of 
probability currents^. We consider a series of configurations a\, o"2, ..cr n , which form a cycle, 
each successive state reachable in a single step from its predecessor (likewise for states a n 
and <Ti): We define the products of the rates around the cycle as: 



il + = C [a 1 G 2 ] C [02 -> CT 3 ] ■ ' • C [<J n -> 0"i] 
Il_ = C [cr 2 -> cri] C [(T3 -> cr 2 ] ■ • • c [(Ti -> cr„] 



(7) 



Detailed balance holds if and only if 



n, = n_ 



(8) 



for all cycles, which is equivalent to saying that the net probability current between any 
two configurations vanishes in the steady state: 
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If the rates violate the detailed balance condition, then there will be non-trivial current 
loops: 
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(10) 



The presence of these current loops is a key characteristic of non-equilibrium steady states, 
and a signal of the microscopic irreversibility of these systems. For a complete and unique 
characterization of a non-equilibrium steady-state one needs to specify both the configura- 
tional probability distribution and the distribution of the probability currents.- The choice 
of the transition rates is not unique. The transition rates we will use (see IHip for our 
two-temperature kinetic model have been previously studied^ and are known to lead to the 
Ising model solution for the equilibrium case. We will examine the specific probability cur- 
rents for our model in order to illustrate the profound differences between equilibrium and 
non-equilibrium systems. 

From a pedagogical point of view, these calculations should provide students with a 
deeper understanding of the profound difference between equilibrium and non-equilibrium 
statistical physics. Starting "from scratch", at the microscopic level, students have the op- 
portunity to learn important statistical physics concepts such as: choice of transition rates, 
balance equations, steady state, probability distributions, equivalence classes and boundary 
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conditions. Students who are familiar with equilibrium statistical physics methods are now 
faced with the challenge of being outside of the traditional framework of the canonical dis- 
tribution, in search of new methods of study. Also, these analytical calculations complement 
straightforward computer simulations of nonequilibrium states. 

In the following section, we introduce our model and present some sample calculations and 
results for an N = 4 lattice. It is quite instructive to see the transformation in configurational 
probabilities as the non-equilibrium condition is "turned on" by gradually allowing two 
reservoirs that permit energy flow in the system to take on different temperatures. 

III. THE TWO TEMPERATURE KINETIC MODEL 

The Ising model was introduced by Lenz in 1925^^, in an attempt to understand the 
nature of phase transitions in ferromagnets. It has become a paradigm of statistical physics, 
and is a common feature of statistical physics classes at the undergraduate level. Building a 
non-equilibrium model related to the Ising model is a good pedagogical approach - students 
already familiar with the Ising model can quickly come to appreciate the novel behaviors 
of systems that are far from equilibrium. Paralleling the standard definition of the Ising 
model, we define our one-dimensional system as follows. Consider a collection of adjacent 
sites numbered i = 1, 2, . . . , N, with site N considered adjacent to site 1, as if the points 
were distributed around a ring (so-called "periodic boundary condition".) Each site i can 
be full or empty, and in order to describe a particular system configuration we define a set 
of occupation numbers rij, each being for an empty site and 1 for an occupied site. We will 
alternately use spin language notation Oi = 2rii — 1, with <j{ = 1 for a full cell and Oi = — 1 
for an empty cell. We assume an even number of sites; periodic boundary conditions imply 
that o"at + i = <j\. The spins are in contact with two heat baths at different temperatures T e 
and T a , in such a way that spins on even lattice sites experience T e and those on odd sites, 
T Q . Imposing T e ^ T Q drives the system out of equilibrium: each heat bath tries to drive the 
system towards equilibrium with the same Hamiltonian but at its own temperature. As a 
result, energy flows from one sublattice to the other and the steady state is a nonequilibrium 
one. 

We endow the spins with nearest-neighbor interactions, according to the usual Ising 
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Hamiltonian: 

H = -JJ2 l a l a l+l , (11) 

where J represents half the energy difference between a state of two adjacent parallel spins 
and one of two adjacent opposite spins. The dynamics is modeled by a generalization of the 
Glauber rates&. The n-th spin is flipped on a time scale of r with a rate 

C n ({O"}) = — (l - y(X n (<7 n+ i + CT n _l)) (12) 

where < 7„ < 1 is related to the local temperature according to 



7n 



W 2J N 

tanh(— — — ) 
k B T e 



n even 

for (13) 



2J 

tanh( — — ) n odd 



' k B T ' 

These rates are invariant under a global spin flip {a} — > {— cr} and under translations by 
an even number of sites a n — > <r n +2j for all n, j = 1,2, ..N (translational invariance modulo 
2). Therefore, the same invariance should hold for the steady state distribution P* ({cr}) , 
namely, P* ({a}) = P* ({-cr}) and P* ({cr n }) = P* ({a n+2j }) This will allow us to reduce 
the number of configurations considered in our model. 

When the two heat baths have the same temperature, T e = T Q = T, the above Hamilto- 
nian and rates define the exactly solvable Glauber model 6 which relaxes to the equilibrium 
state of the Ising model. 

The non-equilibrium version of this one-dimensional kinetic model has been previously 
studied^&^i. The magnetization (average over all spins) and the two-spin correlations (aver- 
age over all pairs of spins) were calculated for both the steady state and the time-dependent 
case.-&i£ So far, there is no compact expression for the steady-state probability distribu- 
tion. Here, we exhibit an exact expression for the full probability distribution, but at the 
price of restricting ourselves to very small systems. 

Below we show the principal steps of the calculations and the most significant results. 

IV. PROBABILITY DISTRIBUTION 

The master equation tells us how a particular configuration evolves in time: 
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N 

d t P({a},t) = £ [cn({cr ln] })P({v [n] },t) -c n ({a})P({a},t)] (14) 

n=l 

where the state {er^} differs from {a} by a flipping of the n-th. spin and the rates are given 
by Eqn. ( 1T21) . We seek the steady state solution: 

P*({a}) = hm^ 00 P({a},t) (15) 

To simplify the problem, we note that configurations which are related by symmetries of 
the dynamics will obviously occur with the same probability. For example, the configurations 

H h+ and + + H — will have the same probability in an N = 4 system, since they result 

from each other by a translation modulo 2. Similarly, H h+ and — I are related by a 

global spin flip. We can therefore define "equivalence classes" , each class consisting of those 
configurations related to one another by a symmetry transformation. We need only solve for 
the probability associated with one member of each equivalence class, thus greatly reducing 
the number of unknowns below the total number of configurations, 2 N . 

Below, we pursue this program for an N = 4 system, present the calculations in some 
detail. We also show the results for an liV = 6 system. The exact stationary probabilities 
will be determined and studied as a function of a single parameter, 7 e for some fixed, suitably 
chosen value of j . Clearly, the equilibrium limit is represented by 7 e = j a . 

For an N = 4 system, there are 6 equivalence classes, numbered (in some arbitrary order) 
i = 0,1,..., 5. The degeneracy di of each class is defined as the number of configurations 
in this class. Pi denotes the stationary probability associated with class i. In Fig. IIVI we 
present the equivalence classes with one representative of each class. For comparison, the 
12 equivalence classes for the N = 6 system are also listed. 

To write the steady state equation for P , we need to identify the "neighbors" of class 
in configuration space, i.e., those configurations which can be reached from Pq via a single 
spin flip, and vice versa: these configurations obviously belong to classes 1 or 2. Similarly, 
class 1 is a neighbor to classes 0, 3, 4 and 5. The master equation leads us to a set of six 
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rate equations: 



2rd t P = 2(l + 7e )P 1 + 2(l + 7 o )P 2 -2(2-7 e - 7o )P 
2rd t P 1 = (1 - 7e )Po + Ps + Pa + (1 + 7e) A - 4Pi 
2r^P 2 = (1 - 7o)Po + P 3 + P 4 + (1 + lo) P 5 ~ 4P 2 
2r<9 t P 3 = 2P a + 2P 2 - 4P 3 
2rd t P A = 2Pi + 2P 2 - 4P 4 

2r^P 5 = (1 - lo )P 2 + (1 - 7e ) Pa - (2 + 7o + 7e ) P 5 . (16) 
Since our probabilities are normalized, we have one additional equation, namely 

5 

l = J2diPi (17) 

The probabilities for the equilibrium case 7 G = 7 e = 7 are given by the Boltzmann factor, 
exp(— j3H). In this case, probabilities differ only if their configurational energies are distinct, 
producing an even greater reduction in the number of distinct cases to consider. With a 
little algebra, one can convert the exponentials into functions of 7 to arrive at: 



Pi = P 2 = P 3 = Pa = P 



0- 



1 + 7) 



P 5 = (18) 

(1+7) 



Using the normalization condition (Eq. [TTj) . we find for Pq the following expression: 

D l(l + 7) 2 



(19) 



u 8(2- 7 2 ) 

These equilibrium results can also be found by solving the rate equations, Eq. El with 
the left-hand sides set to zero, along with Eq. [TTJ 

As we can see, at equilibrium only three different probabilities remain, reflecting the three 
possible values of configurational energy: (i) no broken (i.e., H — ) bonds - class 0; (ii) two 
broken bonds - classes 1, 2, 3, 4 (iii) four broken bonds - class 5. In Fig. Ela we show their 
dependence on 7. All probabilities become equal for 7 — > which corresponds to the infinite 
temperature limit, and all except Po vanish for 7 — > 1, i.e., T —>■ 0. The situation for the 
N = 6 system, shown in Fig. [3la is similar. There are four distinct probabilities, all of which 
are zero at infinite temperature, and three of which vanish at zero temperature. 
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The problem becomes more difficult when equilibrium is not assumed, i.e., 7 Q 7^ 7 e . It 
becomes necessary to solve simultaneously for the probabilities. Steady state solutions are 
achieved by considering the equations, Eq. [2J resulting from the master equation, with the 
time partials all set to zero. The normalization condition, Eq. [T71 implies that the six 
equations of Eq. [2] must be redundant. A solution emerges if one solves the inhomogeneous 
set of linear equations resulting from taking five of the equations from Eq. [2] along with Eq. 
[T71 This is easily done using algebraic software, and produces the results 

8 + 7o 2 -67o7e-3 7e 2 



Pi 

P-2 

Pa 
P 5 



64(2 


- 7o7e) 


8 - 3 7o 2 - 


- 67 7 e + 7 2 


64(2 


- 7o7e) 


8-7o 2 - 


67o7e - 7e 


64(2 


- 7o7e) 


8-7o 2 - 


67o7e - 7e 



64(2 - 7o7e ) 

- 3 7o 2 - 2 7o7e - 3 7e 2 + 8 7o + 8 7e 
64(2 - 7o7e ) 



p 8 + 3 7o 2 + 2 7o7e + 3 7e 2 + 8 7o + 8 7e 

P ° " 64(2 - 7o7 e) • (20) 



We exhibit these probabilities in Fig. [2jb, where we plot each probability vs j e , for 7 D = 0.5, 
letting 7 e run through its full range, to 1. First of all, there are more distinct probabilities 
than for the equilibrium case. The equations show that each equivalence class behaves 
differently, except for the added degeneracy that results from the dynamics we have chosen: 
P3 = P4. The crossings seen at the point where 7 e = 70 simply reproduce the equilibrium 
results for 7 = .5. We also note that Pq remains the most probable configuration. Finally, 
we observe a grouping of curves: The probabilities of configurations which share equal 
configurational energy track each other quite closely, and those associated with different 
configurational energy never cross. One might be tempted to conjecture that this is a generic 
feature of this simple non-equilibrium system. Unfortunately, it does not persist for larger 
system sizes, such as N = 8. This study also allows us to seek configurations, unrelated 
by symmetry, which would nevertheless occur with the same probability. In equilibrium, 
such configurations would all have the same energy, i.e., degeneracies (beyond symmetry) 
are controlled by energy. Far from equilibrium, it is not known which quantity controls such 
degeneracies. The difference between equilibrium and non-equilibrium probabilities is quite 
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dramatic, and any connections between the two are far from obvious. 



V. PROBABILITY CURRENTS 



The presence of current loops marks a fundamental difference between non-equilibrium 
and equilibrium systems, and we here illustrate this for our N = A system. In Section [III we 
defined these probability current loops as: 



J [ a ' -> a] = c 



a 



a 



P*(a 



a 



a 



P*{a) ^ 



(21) 



In the equilibrium case, due to the detailed balance condition, these current loops vanish 
for any pair of configurations in the steady state. For our N = 4 system in steady state, we 
can calculate these currents as: 



Joi = (l-7e)P -(l + 7e 


)Px 


(22) 


J02 = (1 -lo)P - (l + 7o 


)P2 


(23) 


Jl3 = JlA — Pi ~ 


P3 


(24) 


J23 — J24 — P2 — 


P3 


(25) 


Jl5 = (l-7e)A-(l+7 e 


)P 5 


(26) 


J25 = (1 - lo)P2 " (1 + 7o 


)Pk 


(27) 






(28) 


where represents the current between configurations ' 


i" and ' 


j" , where i,j — 0...5. Fig. 


HI shows these currents as a function of j e , for 7 Q = 0.5. 


We can 


see how all currents vanish 



at equilibrium, where 7 e = 7 G = 0.5, in accordance with the detailed balance condition. 



VI. TIME DEPENDENCE 



In principle, the present model can be fully and analytically solved for the time depen- 
dence of the configuration probabilities, given any initial set of probabilities as an initial 
condition. The master equation leads to a set of coupled, first-order differential equations in 
time that can be solved using standard methods. However the size of the system of equations 
grows very quickly with the number of cells, making such a straightforward approach dif- 
ficult and rendering computer simulations a much more reasonable approach. Nonetheless, 
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it is instructive to look at one nontrivial case where complete solutions can be found: the 
N = A case. 

As we have seen, this case allows 2 4 = 16 distinct spin configurations, leading to 16 
distinct time-dependent probability functions. When examining asymptotic time behavior 
the probabilities group into six equivalence classes, leaving a very tractable 6x6 linear 
algebraic system to solve for asymptotic probabilities. For arbitrary times, this simplification 
is not available, since the initial conditions need not obey the symmetries of the asymptotic 
situation. 

It is a straightforward task to set up the 16 time-dependent equations for the probabilities 
that follow from the master equation. They can be simply represented by 

2rd t P = AP, 

where P is a 16 component vector whose components are the configurational probabilities, 
and A is the matrix of coefficients of the probabilities from the right-hand side of the master 
equation. If we assume a time dependence for the vector P of exp (—A) the equation reduces 
to the eigenvalue equation for the matrix A its eigenvalues Oj are proportional to the decay 
coefficients A: A« = aj/(2r). One of the standard algebraic software packages, like Maple, 
can easily identify the eigenvalues for this system, as a function of the parameters 7 e and 
7o, as 0, 2, 4, 6, 8, 2(1 ± y/jey ), 4(1 ± y/j^%) , 6(1 ± y/j^%) ■ The eigenvalues 2 and 6 are both 
twofold degenerate, and 4 is fourfold degenerate. It is worthy of notice that these eigenvalues 
are coincident with the set given by Glauber for the one-temperature case of this model, in 
the limit 7 e = j a . 

The general time-dependent solution is 

P(t) = ^c^exp(-A^), 

i 

where P(t) is a 16 component vector carrying the time dependence of each configuration 
probability, Vi is the eigenvalue of A corresponding to eigenvalue Aj, and the sum runs 
over all eigenvectors. In the case of degenerate eigenvalues, a linearly independent set of 
eigenvectors is chosen. The constants q are fully determined by the initial probabilities: 

C = v- 1 ^ 

— * 

where C is the vector with components q, V is a matrix whose columns are the eigenvectors 

— * 

of A, and Pi is the vector with components equal to the initial probabilities. 
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From the general solution it is clear that the asymptotic time behavior results solely from 
the term corresponding to the A = eigenvalue. The associated eigenvector, normalized, 
coincides with the collection of probabilities derived earlier for the steady-state solution. 

By way of example, we present graphs of time-dependent behavior of the N = 4 system 
in two specific examples. In both cases, the system at t = is fully in the configuration 
In the first, simpler case, we look at the situation in which the two temperature 
baths are at the same temperature, so that the eventual behavior of the system tends to 
equilibrium. Time dependence of the probabilities for this system are exhibited in Fig. 
EJ Because of the simplicity of the initial state and the symmetry introduced by the equal 
temperatures, only five of the sixteen configuration probabilities show distinct time behavior. 
At large values of t, these five probabilities collapse into three values corresponding to the 
equivalence classes exhibited in Fig. IIVI (c). The fact that even in this simple case some 
probabilities grow from zero to a maximum before falling to their asymptotic values indicates 
some nontrivial physical effects are possible. Looking at the same initial condition, but for 
two distinct temperatures (j e = .2,j = .8), we are examining a system that never reaches 
equilibrium, but does eventually achieve a steady state, as discussed earlier. As shown in 
Fig.s [6] and Ui eight of the sixteen configurations evolve distinctly, eventually reaching the 
steady state of the six equivalence classes listed in Fig. HVl a. As for the former example, 
there are configuration probabilities that grow or decay monotonically from their initial to 
final values, but others that initially grow past their steady state values before relaxing 
to them. Since one can describe all physical properties of this system as functions of the 
configuration probabilities, this behavior suggests the possibility of time-dependent peaks 
in some of these properties. Such behavior in this simple system can provide clues in the 
search for interesting physical behaviors in larger systems, which can be investigated using 
computer modeling. 



VII. CONCLUSIONS 



We presented a simple, intuitive way to introduce students to the field of non-equilibrium 
statistical physics using a one dimensional spin system, the two-temperature kinetic Ising 
model. We solved exactly the master equation for a 1x4 system and found the probability 
distribution for both the steady-state and the time dependent case. By comparing the 
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steady-state probability spectrum with the equilibrium counterpart, we were able to see the 
dramatic differences between the two systems. For example, the non-equilibrium steady 
state probability distribution is governed not only by the configurational energy (which 
is the case for equilibrium), but also by other factors, that remain to be explored. The 
time-dependent probability spectrum offers new features that are worth pursuing (e.g local 
maxima that may be indicative of an oscillatory behavior of the system before settling into 
its steady-state). We also emphasized the presence of the probability currents- a defining 
feature of non- equilibrium systems. 

From a pedagogical point of view, this project offers students the opportunity to observe 
and learn about the novel behavior of a driven system, and face the challenges of a lack of 
theoretical framework for far from equilibrium systems. We see this analytical study as a 
starting point for a more complex project. We plan to study bigger system sizes with the 
help of Monte Carlo computer simulations. Besides being a useful and necessary tool in the 
study of various theoretical models, computer simulations are very appealing and accessible 
to undergraduate students. Programming is an important part of any scientific curriculum. 
Research projects based on computer simulations give students an opportunity to bridge 
three core disciplines of their curriculum, physics, computer science and mathematics. Also, 
computer simulations are a intuitive way to introduce students to the non-equilibrium sta- 
tistical physics: using visualization tools, students can "see" how the system goes from an 
equilibrium to a non-equilibrium state by simply switching off and on the extra heat bath. 
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(a) Equivalence classes for 1x4 two temperature 
kinetic model (non-equilibrium). 
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(c) Equivalence classes for 1x4 two temperature 
kinetic model (equilibrium). 



FIG. 1: Equivalence 
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(b) Equivalence classes for 1x6 two temperature 
kinetic model (non-equilibrium). 
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(d) Equivalence classes for 1x6 two temperature 
kinetic model (equilibrium). 



for 1x4 and 1x6 systems 
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(a) Equilibrium 




2J 

FIG. 2: (a) Probability distribution of the 1x4 system as a function of 7 = tanh(- ) in the 

k B T 

equilibrium case, when even and odd sites are in contact with heat baths at the same temperature; 
(b) The probability distribution of the 1x4 system as a function of 7 e in the driven case, when even 
and odd sites are in contact with heat baths at different temperatures, with j = 0.5. 
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2J 

FIG. 3: (a) Probability distribution of the 1x6 system as a function of 7 = tanh(- ) in the 

k B T 

equilibrium case, when even and odd sites are in contact with heat baths at the same temperature; 

(b) The probability distribution of the 1x6 system as a function of m = — in the driven case, 

7o 

when even and odd sites are in contact with heat baths at different temperatures, with 7 G = 0.5. 
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FIG. 4: The probability currents for the 1x4 system as a function of 7 e in the driven case, when 
the even and odd sites are in contact with heat baths at different temperatures, with j D = 0.5. 
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FIG. 5: Time dependence of configuration probabilities in the N = 4 case for 7 e = j D = .5, with 
the ++++ configuration as the initial condition. Vertical axis is absolute probability, time units 
are arbitrary. 
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FIG. 6: Time dependence of configuration probabilities in the N = 4 case for 7 e = .2, 7 D = .8, with 
the ++++ configuration as the initial condition. Vertical axis is absolute probability, time units 
are arbitrary, (a) 
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FIG. 7: Time dependence of configuration probabilities in the N = 4 case for 7 e = .2, 7 D = .8, with 
the ++++ configuration as the initial condition. Vertical axis is absolute probability, time units 
are arbitrary, (b) 
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